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Abstract:  The  sensitivity  of  a  solution  of  a  parametrized  equation 
F(z,  A)  =  0  with  respect  to  the  parameter  vector  A  is  usually  defined 
as  the  change  of  the  state  z  in  dependence  of  X.  In  other  words ,  for 
any  solution  expressible  in  the  form  (?( A).  A)  with  some  smooth  func¬ 
tion  z  —  (z(  A)  the  sensitivity  is  the  derivative  Dz( A).  Typically  the 
solutions  form  a  manifold  M  in  the  product  of  the  state-space  and  the 
parameter  space  and  this  sensitivity  is  available  only  at  those  points  of 
M  where  the  parameters  can  be  used  to  define  a  local  coordinate  sys¬ 
tem.  This  paper  introduces  a  general  sensitivity  concept  which  applies 
at  all  solutions  on  M  and  which  includes  the  earlier  definition.  Some 
general  geometric  interpretations  of  the  new  measure  are  presented  and 
it  is  shown  that  the  sensitivity  analysis  can  be  easily  integrated  into  the 
solution  process.  The  theory  also  suggests  the  introduction  of  a  readily 
computable  second-order  sensitivity  measure  reflecting  the  curvature- 
behavior  of  M.  Two  numerical  examples  illustrate  the  discussion. 
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1.  Introduction 


Mathematical  models  for  many  scientific  and  engineering  problems  have  the  form  of  a 
nonlinear  equation 


(1.1) 


F(z,  A)  =  0, 


1  involving  some  (often  infinite-dimensional)  state  variable  z  and  a  finite-dimensional  vector 

i 

A  of  parameters.  The  computational  tasks  then  include  not  only  the  calculation  of  suitable 
solutions  but  also  the  determination  of  their  sensitivity  under  variations  of  the  parameters. 


In  recent  years  the  literature  on  methods  for  a  sensitivity  analysis  of  specific  problems 

lias  been  growing  rapidly  (see,  e.g.,  [7],  [11]  and  the  references  given  there).  Most  of  this 

*  The  work  was  partly  supported  by  ONR-grant  N-00014-90-J- 1025  and  NSF-grant 
CCR-S907G54.  The  paper  was  begun  during  a  stay  at  tin.  Cenue  for  Mathematical  Analysis 
.n  the  Australian  National  University. 
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work  is  based  on  the  assumption  that  near  a  given  solution  (c0,Ao)  of  (1-1)  the  state  c 
depends  smoothly  on  the  parameter  vector  A;  that  is.  the  solutions  of  (1.1)  can  be  witten 
in  the  form  (z(  A),  A)  with  some  smooth  function  c  =  c(A).  Then  the  derivative  Dz(  A0 ) 
of  ;  at  Ao  is  a  natural  measure  of  the  sensitivity  of  the  solution  (r(Ao).A0).  For  suffiently 
smooth  F  this  sensitivity  measure  has  to  satisfy  the  linearized  equation 

(1.2)  DzF( ~o-  A0 )Dz(  Ao )  +  D\F(zo.  A0 )  =  0. 

If  this  equation  's  explicitly  available  and  D:F(z0,  A0)  is  invertible  then  Dz{  A0)  can  be 
computed  from  (1.2).  Codes  for  solving  (1.2)  for  various  applications  in  molecular  dynamics 
are  cited  in  [7].  When  the  partial  derivatives  DZF  and  D\F  are  not  explicitly  known  we 
may  consider  the  use  of  some  difference  approximations  for  them  in  which  case  there  is  also 
a  need  for  studying  the  influence  of  these  approximations  upon  the  desired  solution  of  (1.2). 
Alternately,  if  the  sensitivity  Dz(\q)/j  in  the  direction  of  a  parameter  vector  fi  is  required 
then  one  may  compute  approximations  of  z{  Ao  +  r/i)  for  several  values  of  the  scalar  r  near 
r  —  0  and  then  apply  some  numerical  differentiation  formula  for  approximating  Dz(Xo)u- 

Generally,  this  definition  of  sensitivity  does  not  reflect  any  of  the  underlying  geometric 
aspects:  moreover,  the  indicated  methods  do  not  take  much  account  of  information  that 
may  be  available  during  the  computation  of  the  original  solution  (co-^o)  itself.  In  fact, 
sensitivity  analysis  is  usually  considered  to  be  a  "post-processing”  technique  which  it 
applied  only  after  a  suitable  solution  has  already  been  found.  The  purpose  of  this  paper  is 
to  give  some  general  geometric  interpretations  of  the  sensitivity  concept  and  to  shew  that 
the  sensitivity  analysis  can  be  easily  integrated  into  the  primary  solution  process  without 
tidding  unduly  to  the  computational  cost.  The  approach  is  based  on  the  application  of  some 
differential-geometric  considerations  which  have  been  shown  earlier  to  be  very  natural  in 
connection  with  parameterized  equations  (1.1).  (see  e.g.  [3],  ISl.  [0]' 


9 


2.  Background 


The  operator  F  in  (1.1)  may  represent  a  system  of  nonlinear  equations  in  several  real 
variables  or  some  boundary  value  problem  in  the  state  variable  c.  In  the  latter  case,  a 
discretization  has  to  be  introduced  for  the  computational  solution  of  the  problem,  and  thus 
in  either  case  we  arrive  at  a  finite-dimensional  equation. 

It  is  useful  to  combine  the  state  variable  c  and  pax'ameter  variable  A  of  (1.1)  into  a 
single  vector  .r  and  hence  to  write  the  equation  in  the  form 

(2.1)  F(x)  =  0. 

Here  we  assume  that  F  :  E  C  Rn  •— » »  Rm .  d  =  n  —  m  >  1.  is  of  class  Cr,  r  >  1  on  an  open 
subset  E  C  Rn  and  0  £  Rm  is  a  regular  value  of  F;  that  is,  rank  DF(x)  =  m  for  all  x  in 
the  inverse  image  F'~1,(0). 

Under  these  conditions  it  is  well-known  that  the  set  of  all  regular  solutions. 

(2.2)  M  =  {  x  6  £;  F(x)  =  0  } 

is  either  empty  or  a  (/-dimensional  Cr-manifold  in  Rn  without  boundary.  We  assume 
always  that  M  ^  0.  For  simplicity  the  tangent  space  TXM  at  any  point  x  of  ill  will  be 
identified  here  with  the  (/-dimensional  affine  space 

(2.3)  TrM  =  {{x,p)  €  {x}  x  Rn\  DF(x)p  =  0}. 

For  numerical  purposes  a  computationally  feasible  scheme  is  required  for  fixing  local 
coordma.t  .w  -  terns  at  a  given  point  xo  €  M.  We  consider  only  orthogonal  coordinate 
systems  and  call  a  linear  map 

(2.4)  V  =  (Vm.V,j)  G  L(Rm  x  /?<*,/?");  VTV  =  Im+d,  VVT  = 
a  local  basis  at  .r()  6  A/  if 

(2.5)  Tx  nker  DF(x0)  =  {0},  T  =  rge  IT  T ^  =  rge  Um. 
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In  terms  of  (2.4),  the  condition  (2.5)  is  equivalent  with 


(2.6)  DF(x0)Vm  €  Isom  (Rm), 
or,  alternately,  with 

(2.7)  {D\?0))  els°m{Rn)- 
We  introduce  the  function 

(2.S)  G  :  Yt(E  -  x0)  C  Rm  *  Rd  ^  Rm,  G(y,  t )  =  F(x0  +  +  Vmy). 

Then,  by  the  implicit  function  theorem  applied  to  the  equation 

(2.9)  G(u,t)  =  0, 

there  are  open  neighborhoods  Sd  C  Rd  of  0  £  Rd  and  S„  C  Rn  of  x0,  respectively,  such 
that  for  any  r  £  Sd  there  exists  exactly  one  solution  y  of  (2.9)  with  xo  +  Vdr  -f  Vmy  £  Sn 
and  that  the  mapping  0  :  Sd  *— *  Rm ,  4>(t)  =  y  is  of  class  Cr  on  Sd .  Evidently,  we  have 
b(0)  =  0  and 

(2.10)  <l>:SdCRd^Rn,  $(T)  =  x0  +  VdT  +  Vm<t>(T),  Vr£Sd, 

is  a  CT-diffeomorphism  from  Sd  onto  M  fl  S„.  In  other  words,  $-1  is  a  chart  of  M  at  xq 
and  we  call  $  the  local  coordinate  map  at  x0  induced  by  the  local  basis  V  and  refer  to  O 
as  the  corresponding  corrector  function. 

For  any  d-dimensional  subspace  T  C  Rn  we  can  choose  an  orthonormal  basis  tq . 

of  Rn  such  that  T  =  span  {?>m+i , . . . ,  vn }  and  then  set  V  -  (?q , . . . .  i'„ ).  Since  the 
condition  (2.5)  depends  only  on  T ,  it  makes  sense  to  speak  of  a  local  coordinate  system 
induced  by  the  subspace  T.  A  point  xq  £  M  is  called  a  foldpoint  with  respect  to  the 
subspace  T  if  the  condition  (2.5)  is  violated;  that  is,  if  T  does  not.  induce  a  local  coordinate 
system. 
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Obviously  (2.5)  always  holds  for  the  d-dimensional  subspace  T  =  her  DF(x o).  In  view 
of  our  definition  (2.3)  of  the  tangent  spaces,  a  local  coordinate  map  induced  by  ker  DF(x0 ) 
shall  be  called  a  tangent-coordinate  map  and  the  corresponding  local  basis  a  tangential 
local  basis.  We  shall  always  denote  tangential  bases  by  U  =  ( Um,Ud )  and  their  induced 
local  coordinate  map  by 

(2.11)  ^  Sd  C  Rd  ^  R\  *(u>)  =  to  +  Udu  +  Vmip(u>),  Vjj  €  Sd, 

Note  that  any  tangential  local  basis  U  is  characterized  by  (2.4)  together  with  the  condition 
DF(xo)Ud  —  0  which  automatically  implies  (2.5). 

As  suggested  by  (1.1),  in  many  applications  the  equation  (2.1)  typically  represents 
some  multi-parameter  problem  which  means  that  some  natural  d-dimensional  parameter 
subspace  A  C  Rn  has  been  identified.  Then  A  defines  a  natural  local  coordinate  system  at 
all  those  points  x  E  M  where  Ax  D  ker  DF(x')  =  {0}. 

3.  First  Order  Sensitivity 

Suppose  that  a  local  basis  map  (2.4/5)  has  been  chosen  at  the  point  x0  E  M  and  that 
(2.10)  denotes  the  induced  local  coordinate  mapping.  Then  the  derivative 

(3.1)  D$( 0)/x  =  Vdfi  +  n  E  Rd , 

represents  the  "change”  of  the  solution  x0  in  the  local  coordinate  direction  /j  E  Rd .  Ac¬ 
cordingly.  it  is  natural  to  define  the  derivative  of  the  corrector  function  <f>  :  Sd  t-+  Rm ;  that 
is.  the  linear  map 

(3.2)  E  e  L(Rd,  Rm ),  T  =  T»d(0). 

as  the  sensitivity  map  at  with  respect  to  the  local  basis  V  and  to  call  the  image  vector 
T//  E  Rm  the  sensitivity  of  the  solution  xc  E  M  in  the  direction  //.  It  is  essential  to  observe 
that  this  definition  of  sensitivity  depends  intrinsically  on  the  choice  of  the  local  basis  V. 


In  applications  where  a  natural,  d-dimensional  parameter  subspace  A  C  Rn  is  given, 
the  sensitivity  map  can  be  defined  at  all  those  points  of  M  where  A  induces  a  local 
coordinate  system.  In  that  case  we  speak  of  the  sensitivity  mapping  at  these  points  with 
respect  to  natural  parameter  changes. 

In  order  to  relate  our  sensitivity  definition  to  that  indicated  in  the  Introduction,  let 
V  =  ( Vm, Vd)  be  a  local  basis  at  x0  with  corresponding  local  coordinate  map  (2.10). 
Recall,  that  y  =  6{t)  is  for  any  r  £  Sd  the  unique  solution  y  of  the  equation  (2.9)  such 
that  Jo  +  E dT  +  I’m 2/  £  S'n •  Evidently,  we  have 

(3.3)  DrG'(0,0)  =  DF(xo)Vd,  DyG(0,0)  =  DF(x0)Vm, 

and  by  differentiation  of  G($(r),r)  =  0  it  follows  that 

(3.4)  DvG(0,0)D  $(0)  +  DrG(0,0)  =  0. 

This  agrees  exactly  with  the  definition  (1.2)  when  applied  to  the  mapping  G  at  a  point 
where  the  natural  parameter  space  A  may  be  used  to  induce  the  local  coordinate  map. 
From  (3. 3). (3. 4),  and  (2.6)  we  obtain  the  explicit  representation 

(3.5)  E  =  -(DF(x0)Vm)-lDF(x0)Vd, 


which  plays  a  central  role  in  many  sensitivity  studies. 

As  a  typical  simple  example  consider  the  cubic 

F  :  R*  <—>  Rl\  F(x)  =  x\  -  j,  j2  -  j3,  Vj  6  R1. 

Since  DF(.t()  )  =  (3.rf  —  j-2.  —  ;r, ,  —  1)  we  see  that  V  =  J3  defines  a  local  coordinate  mapping 
on  M  for  all  j  £  M  with  3.rf  —  j2  ^  0  and  (3.5)  provides  that 


Note  that  jn  =  0  is  a  foldpoint  with  respect  to  V  where  E  indeed  is  not  defined. 
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Suppose  that  l  =  ( Um,Ud )  defines  a  tangential  local  basis  at  xo  with  corresponding 
local  coordinate  map  (2.11).  Then  by  (3.5)  and  the  tangency  property  DF(x0)lJd  —  0  we 
obtain  S  =  0.  In  other  words,  with  respect  to  any  tangential  local  basis  the  sensitivity 
mapping  always  has  the  value  zero. 

Let  \  =  (Vm,Vd)  again  be  any  local  basis  at  3-0  with  corresponding  local  coordinate 
map  (2.10).  Then  the  relation 


(3.6) 


x  —  x0  +  LjT  +  \'m<f>(T)  —  Xo  +  Ud iu  -f-  Unx'4'{^)i 


has  to  hold  for  all  x  in  some  neighborhood  of  x0  on  M.  From  (3.6)  it  follows  that  the 
coordinate  transformation  relating  u;  and  r  is  given  by 


(3.7 


r  =  7(W)  :=  VfUd«  +  VdJ  UmxKu). 


Hence,  from  D4'( 0)  =  0  it  follows  that  Dy(0)  =  VjUd.  Note  that  by  (2.7) 


VI 


DF{ia)\  ( DF(x„)\n  ( DF(x„)DF(za)r  0  . 

~  T  Uj  )  =(  (DF(z0)Vdf  vrUd)ehom{R  . 


and  hence,  because  of  DF{xq)DF{x0)1  <S  Isom  (Rm),  that  D~/{0)  =  VjUd  G  Isom  (Ra 
as  expected. 

From  (3.6)  we  see  that 

d(r)  =  V£Udu  +  Umt(u>) 

which  together  with  Dc( 0)  =  0  provides  the  new  representation 


( 3.S) 


F  =  VjUd(VTUd) 


for  the  sensisitivity  mapping  with  respect  to  V.  The  equivalence  of  (3.5)  and  (3.8)  can 
also  be  established  diiTct.lv.  In  fact,  from  the  relation 


(3.9 ) 


WT  =  vdvj  +  vmv„\  =  /,„ 


-7’ 


i 


it  follows  tli at 


DF(x  0)(VdVdf  +  VmVl)Ud  =  DF(xu)Ud  =  0, 


and  therefore  that 


-(DF(x0)Vm)-xDF(x0)Vd  =  VZ  Ud(VdrUd)~\ 


as  claimed. 


For  the  computation  of  the  sensitivity  map  S  by  means  of  the  representations  (3.5)  or 
(3.8)  we  have  to  solve  an  m  x  m  system  or  a  d  x  d  system  of  linear  equations,  respectively. 
In  practice,  n  is  large  while  d  is  relatively  small.  Thus  the  computational  cost  for  (3.5) 
generally  exceeds  that  of  (3.S).  In  this  connection,  it  is  important  to  observe  also  that 
(3.S)  does  not  depend  on  the  orthogonality  property  Uj  Ud  =  I,i  but  only  on  the  fact  that 
the  columns  of  U(j  span  the  null-space  her  DF{x o).  Indeed,  for  any  Ud  £  L(Rd.  7?")  with 
DF(xo)Ud  =  0  and  rank  I'd  —  d  there  exists  some  cl  £  Isom  Rd  such  that  U  ,i  =  I.  A  A 
whence 


V^UdiVjUd) 


-1 


-  v£  Ud  A  A  ~ 1  ( Vj  Ud)~l  =-  VZUd(VjUd)-\ 


4.  Some  Geometrical  Aspects 

As  before,  at  the  given  point  ,r0  £  M  let  V  =  (V’m,  Yd)  be  any  local  basis  and  U  =  (Um,  I'd) 
a  tangential  local  basis.  Then  by  (3.S)  the  sensitivity  at  xq  in  the  direction  //  £  Rd  with 
respect  to  I’  is  given  by  the  vector 

n  =  Ufi  =  V^Yd{VjUdr\<  £  R". 

Let  v  —  ( I  J  Ud  )  1  //  and  introduce  the  //-dimensional  vectors 

(4.1  )  \>  —  UdV  £  ker  DF(xt) ),  ft  =  Vd/i  £  T  —  rge  Vf/.  b  =  Ymrr  £  Tx. 
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for  which  evidently 


(4.2) 

Then  we  have 

(4.3) 
and 


IHI2  =  MU, 


*  ,,  \\a  ,  t=  cr 


Wjv  =  f‘>  V  =  v 


(4.4 


ft  +  b  =  L  dl  J  v  +  VdVT v  —  i>. 


In  other  words,  v  is  the  unique  vector  in  ker  DF{x0)  for  which  the  orthogonal  projection 
onto  T  equals  fi.  and  a  is  the  orthogonal  projection  of  v  onto  T±  =  rge 

As  a  simple  scalar  measure  of  the  directional  sensitivity  a  it  is  natural  to  compute  the 
norm  ||rr|j2.  For  this  note  that  by  (4.3) 


(4.5; 


and  therefore 


vTv  =  vTuvTv  =  i/ruuTvuTu  =  ||//||22. 


cos  0  = 


-  T  ~ 

f.1 1  V 


IH2IWI2  IMI2 

■’here  0  <  0  <  7t/2  represents  the  angle  between  /)  and  v.  Moreover,  from  (4.4)  it  follows 
that 


w 


(7rz>  =  || cr)|23 ,  <jTft  =  0. 


whence  in  view  of  (4.2) 

( 4.G ) 

and  thus 

(4.7) 


Mh  =  IMlh  +  IHI/. 


e  2  =  -  2n 
cos  z0 


1  5  -  iil/2ii/<ii2  =  si“ e 


[1  -  sin  20] 


1/2  He  II 2 
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Evidently,  for  0  —  rr/2  the  condition  (2.5)  is  violated  and  obviously  this  occurs  exactly 
when  ./■()  is  a  foldpoint  with  respect  to  T.  Thus  the  sensitivity  provides  a  measure  of  the 
suitability  ot  the  local  coordinate  system  induced  by  T  for  representing  .1/  locally  near  the 
current  point  .r(). 

In  the  special  case  d  —  1  of  one-dimensional  manifolds  the  angle  0  is  explicitly  known. 
In  fact,  for  d  =  1  suppose'  that  DF(sq)u  =  0,  |juj|9  =  1  and  that  T  =  O', „,»■).  «’  £ 
peM  =  1.  is  a  local  basis  at  X(>.  Then  v,  //  £  Rx  and  v  and  /}  are  n- vectors  in  the  direction 
of  n  and  ir.  respectively;  that  is.  we  have 

T 

cos  0  =  <r  a. 


Generally,  the  distance  between  any  two.  equi-dimensional  linear  subspaces  Si  and  S< 
of  R"  defined  by 

dist  (  S,  ,S>)  —  || P)  -  P2\\2 

wh'To  P,  is  the  orthogonal  j>rojection  onto  S, .  i  --  1.2.  In  the  case  of  one-dimensional 
spaces  it  is  easily  seen  (see  e.g.  [4])  that 

dist  ( S i .  S_> )  =  sin  0 


where  0  is  the  angle  between  the  two  subspaces.  Thus  for  d  =  1  and  with  To  =  ker  DF(.i'o) 
we  mav  write  l  4.7 )  as 


„  „  <list  (T.Ta) 

(4.X  hr  ,  =  - ~ -  t/  . 

‘  [1  -  dist  (T.To)2]1/2 

Moreover,  since  in  this  case  E  £  L{  /?'./?'"  ),  m  =  n  —  1.  and  //£/?’.  (4.S)  implies  that 

dist  ( T,  To  ) 

U.fh  EL  =  - - - . 

-  [1  -dist  (T.Toryr* 

It  turns  out  that  (4.9)  holds  generally  for  manifolds  of  any  dimension  d.  In  order  to 
soe  this  let 


M.10)  .4  1  ( \  '/  l ',( )13  =  diag  (cos  0] . ros  0j).  0  <  0j  <  .  .  .  0,/  <  tt/2. 
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be  the  singular  value-  decomposition  cf  Vj  Uj  and  set  again  To  =  her  DF{tq).  Then  it  is 
well  known  that 


f4.il)  dist  ( T,  T0 )  =  sin  , 

(see  e.g.  [GYLj.  p24).  By  (4.G)  we  have 

s^iikHiivyoa-viik-MiT 


and  hence 


=  lllaX  |//||._JIKl'/t:«/)  V  11-2 2 


- 1  = 


:rf  C'rf)"  ll2” _  l- 


Bv  (4.10)  we  sor  that 


and  thus 


ll(,'<'Tt'rl|l*  =  JSre; 


ivii  2  _  _ I _ 1  _  sin 


1 J  cos  20(i  1  —  sin  20rf 

which  together  with  (4.11)  gives  (4.9),  ar  claimed. 


5.  Second  Order  Sensitivity 

Let  V  —  ;  £',/)  be  a  tangential  local  basis  at  6  At  and  y  the  corresponding  corrector 

function.  Then  we  saw  that  the  sensitivity  map  S  =  Dir'(0)  is  zero,  and  lienee  it  is  of 
interest  to  consider  s  lie  second  derivative  D2v{ 0).  For  this  we  assume  from  now  on  that 
our  mapping  F  is  at  least  of  class  Cr  with  r  >  2. 

The  natural  inner  product  of  R"  induces  a  Riemannian  structure  on  the  manifold 
M.  In  [G]  an  algorithm  was  developed  for  computing  the  second  fundamental  tensor  on 
A/.  Moreover,  it  was  shown  that  for  suitable  bases  of  the  tangent  and  normal  spaces 
this  tensor  equals  the  desired  second  derivative  of  the  corrector  function  for  tangential 


11 


local  coordinate  systems.  We  shall  not  detail  this  approach  here  but  sketch  only  the  basic 
numerical  method  of  [6]  as  it  applies  to  the  computation  of  the  corrector  function. 

Let  V  —  ( I’m '  L'd )  be  any  tangential  local  basis  at  x0  €  M  with  (2.11)  as  the  induced 
local  coordinate  map.  For  any  given  //  G  Rd  with  [|/i||2  =  1  we  have  T(//i)  G  Sd  for  all  real 
t  in  some  interval  J  =  {—(.().  Then 


(5.1)  £  :  J  M,  £(  t )  —  Xo  +  tUdfi  +  L mil'itfj),  t  G  J, 


defines  a  path  on  M  through  x0  which  has  at  x0  the  tangent  vector  {x0}  x  £'(0)  = 
{.ro}  x  C>.  Moreover,  we  have  £"( t )  =  /j.)  and  thus  £"(0)r£'(0)  =  0. 

Let  i]  :  Jo  >— >  .7  be  the  transformation  of  the  arclength  >  of  the  path  to  the  parameter 
/.  and  define  C  :  .Ju  >->  M  by  C(.s)  =  (Jr)(s)),  s  G  Jo-  Then  ||c)‘,(0)||2  =  1  ensures  that 
?/'(())  =  1  and  from 


0  =  2C"(MTC'(M  =  VUV/(s)|in^))l|22  +2r/(s)2^"(7?(s))T^(^(5)) 


it  follows  that  //"(0)  =  0.  Therefore  we  have  ("(0)  =  £”(0)  and  hence 


h  =  llf"(0)|l2  =  ||Cm/JV(0)((l,/i)||2 

is  the  curvature  of  the  path  if  at  x0  and  the  direction  of  £"(0)  equals  that  of  the  principal 
normal  of  £. 

With  some  sufficiently  small  1 0  G  J  set  u>\  =  =  —  to[i  and 


o.2) 


J'i  —  -To  T  L  d^'t  +  U rn  t/’(ce,),  i  —  1,2, 


If  the  triangle  defined  by  the  points  x0,x1,X2  is  non-degenerate  then  the  curvature  of  the 
circumscribing  circle  is  given  by  Heron’s  formula 

k-  =  -y- [.s ( ,s  -o)(.s  -  b)(s  -  c)]1/2 
anc 
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where  s  =  (a  +  b  +  c)/2,  a  =  ||xi  —  xo||2,  b  =  ||x2  —  Xo||2.  With  the  normalization  ac  =  a/c , 
bc  =  h/c  this  can  be  re-written  as 

(5.3a)  k  =  -[—  +  — ](1  -  <5)^2sin  a, 

C  CLc  CLc 

where 

(5.3 6)  8  =  ac  —  bc,  a  —  arcos  7,  7  =  - — . 

ac  +  oc 

When  1  —  *  falls  below  the  machine  precision  then,  in  floating  point  arithmetic,  a  will  be 
zero  and  we  set  k  =  0.  Clearly,  if  to  tends  to  zero  then  the  circle  through  x0,xi,x2  tends 
to  the  osculating  circle  of  the  path  at  xo  and  hence  k  becomes  the  curvature  ko-  Thus  our 
algorithm  produces  an  approximation  k  of  k0. 

If  k  is  not  zero  then  an  approximation  fi  of  the  principal  normal  of  the  path  at  x0  can 
be  generated  by  orthogonalizing  the  sum-vector  (i\  —  x0)  -f  (z2  —  xo )  with  respect  to  UjH 
and  then  normalizing  the  result.  In  other  words,  we  apply  the  algorithm 

(5.4)  ft  :=  (xi  -  x0)  +  (x2  -  x0);  h  :=  n  -  [{Udn)Trt]Udn\  n  :=  n/||n||2. 

Thus  altogether 

(5.5)  UmD2t{0)(n,  f*)  ~  kh. 

The  bilinearity  and  symmetry  of  UmD2xk( 0)(/x,  /*)  implies  that  for  any  vectors  /j  1 .  //2  £ 
R'1  and  with  u  =  // j  +  p2  we  have 

UmD2 ir(0)(fji.  //2)  =  ~[UmD24iO)(i/,i/)  -  UmD2xl>(0)(ni,  fix )  —  Z7mZ>2?/’(0)(//2, /i2 )] 

Hence  by  applying  the  above  algorithm  three  times  we  can  compute  also  £’„,D2t/’(0)(//i ,  //2 ). 
Generally,  for  the  evaluation  of  this  derivative  for  arbitrary  pairs  of  vectors  //i,//2  €  Rd 
we  need  to  compute  only  the  d(d  +  l)/2  terms  UmD24,(Q)((J-n  dj ),  1  <  ?  <  j  <  d  for  some 
basis  , . . .  ,fid  of  Rd ,  (see  [6]). 
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Now  let  V  =  (l’m,  I’d )  be  any  other  local  basis  at  r0  E  M  with  corresponding  corrector 
function  0.  Then  it  follows  from  (3.7)  and  D~f(0)  =  VjUd  that  for  any  /ii,//2  E  Rd 

VmD24>mvJ'UdliUVfUdli2)  +  [Vd  +  VmE]D27(0)(/i1,/£2)  =  UmD2il>(  0)(m,//2). 

Hence  with 

D2n(Q)(Hul*‘>)  =  VlUmD24’(0)(nUfi2) 
and  (3.9)  it  follows  readily  that 

(5.G)  \ D2 d( 0 )( Vj Udi-i i ,  Vj Udn2 )  =  [VmV£  -  VfnEVrfT]C7mD2^(0)(/i1,/i2). 

Xow  note  that  by  (3.8) 

(5.7)  Vmv£  -  VmZVf  =  Vm\r£P,  P  =  In~  Vd{VdTUd)-'V? , 

Since  PYm  =  I’m  and  PUd  =  0,  we  see  that  P  E  L(Rn,Rn)  is  the  projection  onto  rge 
parallel  to  rge  Ud.  Moreover,  it  follows  that  VmV£P  =  (/n  —  VjVd)P  —  P  whence  by 
(5.6) 

(5.S)  VmD24>mvIUdiiUVjUdii2)  =  PUmD2^( 0)(/i,,/i2),  /',  €  Prf,  *  =  1,2, 
where  the  right  hand  side  can  be  calculated  by  applying  the  above  algorithm. 

6.  Some  Example  Problems 

A  simple  0-dimensional  model  of  the  global  energy  balance  of  the  earth’s  climate  has  been 
discussed  in  [2]  and  [5]  (sec  also  [10]).  In  dimensionless  form  it  can  be  written  as 

(0.1)  =  ji{l  -  a{r))  -  ear4 ,  a(r)  =  max(min(f>  -  pr, amin), amax). 

Here  the  state  variable  r  is  a  dimensionless  temperature  defined  as  the  quotient  T/T0  of 
an  average  surface  temperature  T  of  a  spherical  planet  and  a  reference  temperature  To- 
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The  function  a  is  a  model  of  the  planetary  albedo  assumed  to  depend  principally  on  the 
extent  of  the  surface  snow  cover.  All  constants  are  dimensionless  and,  in  particular,  c 
characterizes  the  global  heat  capacity.  The  parameter  e  defines  the  effective  emissivity 
which  incorporates  the  effect  of  water  vapor,  carbon  dioxide,  dust,  etc.,  upon  terrestrial 
radiation  while  fi  characterizes  the  variation  of  the  solar  radiation  from  the  reference  value 
f.i  =  l.  The  model  assumes  that  all  quantities  are  annually  averaged  to  remove  the  seasonal 
cycle. 

We  use  the  reference  temperature  7o  =  288. 6°A’  of  [3]  and  the  constants  o  =  1.14. 
6  =  2.8.  Ckm =  0.2,  amaI  =  0.8.  As  in  [5],  for  the  computation  the  albedo-definition  is 
replaced  by  a(r)  =  3(b  —  pr)  where 

0.2  if  s  <0.2-$, 

0.2 +  (s- 0.2 +  6)2/46  if  \s  —  0.2|  <  6, 

,!l(s)=  7  if  0.2  +  6  <  s  <  0.8-6. 

s  -  (s  -  0.8  +  6)746  if  |s  —  0.8|  <  6, 

.  0.8  if  s  >  0.8  +  6. 

with  6  =  10_  i.  Then  the  equilibrium  configurations  of  the  model  (G.l)  are  the  solutions 
of  the  equation 

(0.2)  f(x)  =  -aeT4+(i[l-fi(b-pT)\=0,  Vi  =  (r,e,;;.p)7’6f?4. 


and  these  form  a  three-dimensional  manifold  M  in  R4 . 

The  equilibrium  configuration  x  €  M  is  (asymptotically)  stable  if  t](t)  <  0  where 


df3  , 

T]{x)  =  (I— p  -  4 aer  , 
as 


is  the  derivative  of  /  with  respect  to  r.  Moreover, the  foldpoints  on  M  with  respect  to  the 
space  of  the  three  parameters,  are  the  points  where  rj(x)  vanishes.  These  foldpoints  form 
the  two-dimensional  sub-manifold 


A/o  —  { x  €  A4 ;  x  —  ( 


4(6-1) 


27  pp4 

256a(6  —  l)3 


,  ft,  p)T  } 
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of  M .  Observe  that  on  Mq  the  albedo  has  the  constant  value  0.4.  Table  G.l  lists  some 


points  of  Mq- 


T 

P 

e 

1.4 

1.714 

0.1370// 

1.3 

1.846 

0.1843 p 

1.2 

2.000 

0.2538/i 

1.1 

2.182 

0.3595// 

1.0 

2.400 

0.5263// 

0.98 

2.449 

0.5706// 

0.96 

2.500 

0.6197// 

0.94 

2,553 

0.6741/z 

0.92 

2.609 

0.7347 // 

0.90 

2.666 

0.8022// 

0.88 

2.727 

0.8776 // 

0.86 

2.791 

0.9622// 

Table  6.1  Foldpoint  Submanifold 


The  current  climate  ,t0  €  M  is  assumed  to  have  the  components  r0  =  0.9893, //0  = 
1,  e0  =  0.7071,  po  =  2.6  where  r0  corresponds  to  an  absolute  temperature  of  285. 5°  A'  and  po 
to  the  dimensioned  slope  0.009  for  the  albedo  decrease  used  in  [10].  Since  i](x o)  =  —0.5220 
this  current  point  is  stable. 

In  the  cited  articles  interest  has  focussed  on  the  path  on  M  through  .To  toward  Mq 
when  //  is  allowed  to  decrease  while  e  =  eo  and  p  =  po  are  held  constant.  Some  points  on 
this  path  and  the  corresponding  values  of  the  first-order  sensitivity  a  =  (a j,  <72,  <t3)  are 
given  in  Table  6.2. 


P 

T 

V 

<72 

03 

1  a\  7 

mam 

0.9894 

-0.5226 

-2.090 

1.478 

-1.893 

3.184 

mss 

0.9816 

-0.4500 

-2.352 

1.672 

-2.182 

3.618 

0.99 

0.9730 

-0.3703 

-2.760 

1.971 

-2.628 

4.290 

0.985 

0.9623 

-0.2733 

-3.578 

2.569 

-3.522  . 

5.639 

0.98 

0.9501 

-0.1658 

-5.605 

4.044 

-5.732 

8.979 

0.9775 

0.9411 

-0.8738(-l) 

-10.23 

7.402 

-10.770 

16.60 

0.9758 

0.9307 

0.2462(-3) 

3474. 

-2518. 

3780. 

5718. 

Table  6.2  Decreasing  Solar  Radiation 
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These  data  show  clearly  how  the  sensitivity  increases  as  the  points  approach  the  foldset 
Mq.  Since  e  incorporates  the  so-called  greenhouse  effect  it  is  also  interesting  to  see  the 
effect  of  decreasing  e.  Table  6.3  lists  some  points  on  the  path  on  M  through  Xq  when  e  is 
allowed  to  decrease  while  p  =  and  p  =  po  are  held  constant. 


e 

T 

V 

o\ 

<r2 

<73 

m2 

0.7 

1.003 

-0.6236 

-1.852 

1.297 

-1.609 

2.775 

0.68 

1.037 

-0.8549 

-1.540 

1.047 

-1.213 

2.223 

0.66 

1.066 

-1.046 

-1.408 

0.9290 

-1.019 

1.971 

0.64 

1.093 

-1.215 

-1.341 

0.8582 

-0.8998 

1.829 

0.62 

1.120 

-1.371 

-1.308 

0.8110 

-0.8170 

1.743 

0.60 

1.146 

-1.517 

-1.296 

0.7776 

-0.7555 

1.690 

0.58 

1.172 

-1.656 

-1.298 

0.7530 

-0.7078 

1.659 

Table  6.3  Decreasing  Emissivity 


As  a  second  example  we  consider  the  following  model  of  an  exothermic  chemical 
reaction  with  convective  transport 


-u"  -f  vu'  =  (1  —  u)exp(ct 


A 


),  0  <  t  <  1,  u(0)  =  u(l)  =  0,  a  =  12 In  10, 


(1  +u) 

discussed  in  [1].  We  use  the  upwind  discretization 

-(1  +  uh)zi- 1  +  (2  +  vh)z{  -  zl+i  =  h2{  1  -  Zi)exp(o  - 
l  —  1,  .  .  .  ,  7T,  ^0  1  0. 


A 


(1  +  zi) 


), 


(6.3) 


on  a  uniform  mesh  with  step-size  h  =  l/(k  +  1).  For  k  =  9  the  2-dimensional  solution 
manifold  of  (6.3)  was  triangulated  by  the  algorithm  of  [9]  in  a  neighborhood  of  the  point 
./'o  =  (*0'A.i/)  €  1  for  which  A  =  23.907,  u  =  999.978  and  the  state  vector  2  = 

(z, . 9)7  has  the  components 


2,  =  4.6067(— 3), 
z4  =  2.2570(— 2), 
27  -  5.4582(— 2), 


22  =  9.7851(— 3), 
25  =  3.0790(— 2), 
28  =  7.5449(— 2), 


23  =  1.5694(— 2), 
26  =4.1012(— 2), 
29  =  1.3374(  — 1), 
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For  the  first  23  computed  nodal  points  of  the  triangulation  Table  C.4  lists  the  values 
of  the  last  (and  largest)  component  29  of  the  state  vector  and  the  two  (natural)  parameters 
v  and  A  together  with  the  Euclidean  norm  of  the  sensitivity  E  at  each  point  with  respect 
to  the  parameters. 


node 

X9 

u 

A 

IIS  1, 

UJ) 

1.70238(-1) 

23.9573 

999.957 

0.398724 

(2,7) 

1.70246(-1) 

23.9573 

1000.00 

0.398640 

(1,6) 

1.67282(-1) 

23.9494 

999.935 

0.435477 

(2,6) 

1.55214(-1) 

23.9248 

999.978 

0.680829 

(3.6) 

1.67276(-1) 

23.9493 

1000.02 

0.435561 

(1,5) 

1.50738(-1) 

23.9187 

999.913 

0.851227 

(2,5) 

1.50739(-1) 

23.9187 

999.957 

0.851212 

(3,5) 

1.50739(-1) 

23.9186 

1000.00 

0.S5119S 

(4,5) 

1.50740(-1 ) 

23.9186 

1000.04 

0.851183 

(1,4) 

1.30793(-1) 

23.9070 

999.892 

16.1791 

(2,4) 

1.33738(-1) 

23.9072 

999.935 

8.80483 

(3,4) 

1.33741(-1) 

23.9072 

999.978 

8.79455 

(4,4) 

1.33743(-1) 

23.9072 

1000.02 

8.78402 

(5.4) 

1.30793(-1) 

23.9069 

1000.07 

16.1566 

(2.3) 

1.15878(-1) 

23.9132 

999.913 

1.07936 

(3,3) 

1.15878(-1) 

23.9131 

999.957 

1.07939 

(4,3) 

1.15879(-1) 

23.9131 

1000.00 

1.07942 

(5,3) 

1.15880(-1) 

23.9130 

1000.04 

1.07945 

(3,2) 

1.13146(-1) 

23.9157 

999.935 

0.922609 

(4,2) 

1.13132(-1) 

23.9156 

999.978 

0.921912' 

(5,2) 

1.13119C-1) 

23.9156 

1000.02 

0.921225 

(6,1) 

9.43086(-2) 

23.9451 

999.957 

0.448821 

(7,1) 

9.43010(-2) 

23.9451 

1000.00 

0.448722 

Table  6.4  Sensitivities  at  Triangulation  Nodes 


The  nodes  are  indexed  by  a  pair  of  integers  (i,j)  where  (3,4)  is  the  above  given  point 
./•()•  They  form  30  triangles  specified  by  the  following  node-triples 

{(i-j)AiJ  +  1),(?:  +  ((*,>),(*  +  hj  -  1),(*  +  1 J)), 

1  <  i,  i  +  1  <  4,  1  <  j,j  +  1  <  7,  5  <  i  +  j  <  8. 

In  Table  6.4  the  nodes  (i,4),  i  =  carry  a  significantly  larger  value  of  ||E||2. 

This  indicates  that  this  line  of  nodes  should  be  near  a  foldline  with  respect  to  the  natural 
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parameter  space.  This  is  indeed  the  case.  A  foldpoint  computation  with  a  standard 
augmentation  of  the  mapping  (6.3)  started  from  x0  converged  very  quickly  and  from  the 
resulting  point  a  continuation  process  readily  produced  the  expected  foldline.  Table  6.5 
gives  the  values  of  zg,  A,  and  v  of  some  of  the  computed  points  on  that  foldline. 


Xg 

V 

A 

1.30771(-1) 

23.9074 

999.578 

1.30773(-1) 

23.9071 

999.848 

1.30773(-1) 

23.9070 

999.938 

1.30774(-1) 

23.9070 

999.978 

1.30774(-1) 

23.9069 

999.988 

1.30774(-1) 

23.9069 

1000.02 

1.30775(-1) 

23.9068 

1000.11 

1.30777(-1) 

23.9065 

1000.38 

1.30780(-1) 

23.9061 

1000.76 

Table  6.5  Foldline 

These  results  certainly  show  that  the  sensitivity  norms  represent  an  excellent  indicator 
for  the  existence  of  nearby  foldpoints  which  could  be  used  effectively  in  any  computational 
study  of  the  characteristic  properties  of  the  solution  manifold. 
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